/* Original version: Tim Simcoe 2007-02-20 */
/* Updated: Miikka Rokkanen 2012-08-10 */

/*   This program caculates a "Robust" Covariance Matrix for the 
Poisson QML model with conditional fixed effects. Formulas were
obtained from Wooldridge (1999), Journal of Econometrics  */

capture: program drop xtqmlp
program define xtqmlp, eclass byable(onecall) sort
version 12.0

syntax varlist(fv) [if] [in], [I(varname num) FE IRr CLuster(varname)]
marksample touse

** Make sure user specifies a Group Variable (j)
if length("`j'") == 0 {	
	local j = "`_dta[iis]'"
	if length("`j'") == 0 {
		di in red "you must specify a group variable for the panel"
		exit 198
	}
}

** Check to see that Groups (j) are nested within Clusters

if length("`cluster'") > 0 { 
	quietly {
		bysort `j' : gen flag = `cluster'!=`cluster'[_n-1] & _n>1
		sum flag
		drop flag


	}
	if r(mean)>0 {
		di in red "group variable (j) must be nested within clusters"
		exit 198
	}

	quietly {
		regress `varlist' if `touse', cluster(`cluster')
		display " regress `varlist' if `touse', cluster(`cluster')  `e(N_clust)'"	
		local cs =`e(N_clust)'
	}

}


if "`irr'" == "irr" {
xtpoisson `varlist' if `touse', fe i(`j') irr	
}
else {
xtpoisson `varlist' if `touse', fe i(`j')
}

/* quietly{
display "xtpoisson `varlist' if `touse', fe i(`j')"
}
*/

display "Calculating Robust Standard Errors..."
quietly {

	tempvar mu mubar ybar score
	
	matrix b = e(b)
	matrix V = e(V)			
		
	qui predict double `mu' if e(sample), xb
	qui replace `mu' = exp(`mu') if e(sample)
	qui egen double `mubar' = mean(`mu') if e(sample), by(`j')
	qui egen double `ybar' = mean(`e(depvar)') if e(sample), by(`j')
	qui gen double `score' = `e(depvar)' - `mu'*`ybar'/`mubar' if e(sample)
	
	if length("`cluster'") > 0 {
		_robust `score' if e(sample), variance(V) cluster(`cluster') minus(0)
	}
	if length("`cluster'") == 0 {
		_robust `score' if e(sample), variance(V) cluster(`j') minus(0)
	}

	/* Calculate Wald Chi2 Statistic (first get rid of omitted variables) */

	_ms_omit_info b
	local cols = colsof(b)
	matrix noomit =  J(1,`cols',1) - r(omit)

	mata: newV = select(st_matrix("V"),(st_matrix("noomit")))
	mata: newV = select(newV, (st_matrix("noomit"))')
	mata: st_matrix("newV",newV)
		
	mata: newB = select(st_matrix("b"),(st_matrix("noomit")))
	mata: st_matrix("newB",newB)
	
	local dof = colsof(newB)
	
	mat chi = newB * inv(newV) * newB'
	local chi2 = trace(chi)
	local pval = chi2tail(`dof',`chi2')	
	
	ereturn local cmd ="xtqmlp"
	ereturn local predict ="xtqmlp_refe_p"
	ereturn scalar chi2 = `chi2'
	ereturn scalar p = `pval'										   
	if length("`cluster'") > 0 {
	ereturn scalar N_clust =`cs'
	}
	ereturn repost b = b V = V 
	
}

if "`irr'" == "irr" {
	ereturn display, eform(IRR)
}

if "`irr'" != "irr" {
	ereturn display
}

di in green "Wald chi2(" in yellow `dof' in green ") = " in yellow %8.2f `chi2' /*
*/ "                                " in green "Prob > chi2 = " in yellow %8.4f `pval'
if length("`cluster'") > 0 { 
di in green "Number of Clusters: " in yellow `cs'
}

end



